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Abstract 

Slow dynamics in a fluid are studied in one of the most basic systems possible: polydisperse hard spheres. 
Monodisperse hard spheres cannot be studied as the slow down in dynamics as the density is increased is preempted 
by crystallisation. As the dynamics slow they become more heterogeneous, the spread in the distances traveled by 
different particles in the same time increases. However, the dynamics appears to be less heterogeneous than in hard- 
sphere-like colloids at the same volume fraction. The particles which move least far in a characteristic relaxation 
time and, particularly, the particles which move farthest in the same time are clustered, not randomly distributed 
throughout the sample. We study the dynamics at three different widths of the distribution of diameters of the hard 
spheres. For each width, the relaxation time is the same function of the compressibility factor, suggesting that this 
determines the relaxation time for hard spheres. 



1 Introduction 

The study of hard spheres has taught us much about the 
equilibrium and dynamic behaviour of dense fluids and 
crystals. We have learnt that crystallisation can be a re- 
sult of the more efficient packing in the crystal of atoms or 
colloids, and that diffusion is inhibited in dense fluids by 
repeated collisions of particles with their neighbours. One 
outstanding problem is that of understanding and perhaps 
trying to predict, the dramatic slowdown of the dynamics 
of fluids and the formation of a glass; see Refs. 1-3 for 
introductions to slow dynamics and glasses. An obvious 
thing to do is to study this phenomenon in hard spheres. 
This has been done and the results are presented here. 

Unfortunately, monodisperse hard spheres crystallise 
readily, preventing the study of hard spheres at densities 
high enough to observe very slow dynamics. Crystallisation 
is inhibited if the spheres are not all the same size. Speedy 
studied a binary mixture of hard spheres of two diame- 
ters, 4 which did not crystallise. Here mixtures of spheres 
with continuous ranges of diameters are studied, these are 
called polydisperse mixtures. This will enable us to com- 
pare with experiments on hard-sphere colloids, which are 
inevitably on polydisperse particles,. 5-7 Monte Carlo sim- 
ulations of polydisperse hard spheres have been performed 
by Doliwa and Heuer. 8,9 

In this contribution we will quantify the slow down in 
the dynamics, characterise these slow dynamics, and assess 
the effect of having spheres of different diameters present. 
The slow down is measured by calculating a relaxation 
time. The dynamics are characterised by calculating their 
deviation from what we would expect if each particle was 



diffusing independently. When the deviations are large we 
call the dynamics heterogeneous. When the dynamics are 
heterogeneous, the distribution of speeds and relaxation 
times of individual particles is much broader than the Gaus- 
sian function obtained when particles are diffusing indepen- 
dently. There are many more particles which have traveled 
much farther than the mean than we would expect for in- 
dependent diffusion of each particle. These fast-moving 
particles are not distributed randomly in space, they are 
concentrated in clusters. In other words, fast particles tend 
to be surrounded by other fast particles. 10-13 

Three different widths of the distribution of diameters 
are studied. The polydispersity width is characterised by 
the standard deviation of the distribution in diameters. For 
the three distributions this is 9, 14 and 20%. As in ex- 
periments on hard-sphere colloids 14, 15 if the distribution 
is narrower the spheres can crystallise. Surprisingly, al- 
though increasing the polydispersity at fixed volume frac- 
tion speeds up relaxation, the relaxation time is, within 
simulation resolution, independent of the width of the dis- 
tribution at fixed compressibility factor. The relaxation 
time appears to depend only on the compressibility factor. 
We also study the effect of the heterogeneity in diameters 
on the heterogeneity in the dynamics. 

In the next section we briefly describe our, standard, 
methodology. Section 3 outlines when our simulations crys- 
tallised and quantifies how much polydispersity in the di- 
ameters is required to avoid crystallisation. Section 4 dis- 
cusses a thermodynamic property, the pressure. The last 
but one section contains our results for the dynamics, and 
the last section is a discussion. All results are for fluids at 
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equilibrium not glasses. The focus is on the relaxation of 
the structure of the fluids, diffusion coefficients and other 
properties of motion at times large with respect to the re- 
laxation time are not calculated. 

2 Simulation methodology 

Conventional NpT Monte Carlo and NVE molecular dy- 
namics techniques are used. 16 N, p, V, E, and T are the 
number of particles, the pressure, the volume, the energy 
and the temperature, respectively. The only nonstandard 
part concerns the polydispersity of the system, the presence 
of spheres of different diameters. In theoretical treatments 
of polydispersity first the thermodynamic limit N — > oo, is 
taken and then the number of components present is taken 
to oo. This gives rise to a continuous distribution of sizes of 
spheres 17 ' 18 with a number density px(a)da of spheres in 
the range a to er + der. p = N/V is the total number density 
of spheres. The width of the distribution of diameters can 
be characterised by a (dimcnsionless) width parameter w. 
The larger w is the broader the distribution of sizes present 
in the mixture. In the limit w^Owe recover a monodis- 
perse system. A simple functional form for the distribution 
function x is the hat function: 

C a < a(l - w/2) 

x(a) = I (aw)- 1 er(l - w/2) < a < a(l + w/2) , 

[ a >a(l + w/2) 

(1) 

where a is the mean diameter. The standard deviation 
divided by its mean, of the distribution x(a), s — w/\/Y2. 

Simulations are of course of a finite system — we cannot 
take the thermodynamic limit in a simulation. For finite 
N, x(a) cannot be a continuous function. For our simu- 
lations we generate a set of N diameters by generating a 
set of N pseudo-random numbers in the range a(l — w/2) 
to <f(l + w/2). Ideally we would average over many differ- 
ent realisations of the polydispersity by performing com- 
puter simulations for many different sets of the N diame- 
ters. However, as our systems have long relaxation times 
and so require long runs this is not feasible. Thus all the 
results we will present for a given polydispersity width will 
be for a single system of 1372 particles with one set of the 
1372 diameters. We have performed some simulations with 
different sets of the diameters but the same number of par- 
ticles and the results did not change significantly. We have 
also performed simulations with 256 and 500 particles and 
obtained results which did not differ significantly with the 
exception of some results for the dynamics. We discuss 
these differences below. 

NpT Monte Carlo simulations are used to start with 
and to compress the fluid phase up to the required den- 
sity. NVE molecular dynamics are used to equilibrate at 
a particular density and to obtain all the averages shown, 
including the pressure. All simulation results are for fluids 
at (possibly metastable) equilibrium, not for glasses. Our 
results are in reduced units. We use the mean diameter 



a as our length scale. The reduced pressure p r is then 
p r = Zpa 3 , where p is the number density of particles and 
Z = pV/(NkT) is the compressibility factor. The results 
are all averages over 8 runs. We use the variation between 
the results of the individual runs to estimate statistical er- 
rors. In no case are they significantly larger than the plot 
symbols used in the figures. 

A time scale can obtained by the ratio of a distance to 
a velocity. The velocity we choose is the root mean square 
velocity along one of the 3 Cartesian axes, of a particle 
with a diameter equal to the mean diameter. This veloc- 
ity is equal to (fcT/m) 1 / 2 , where m is the mass of particle 
with diameter a. The time scale is thus aijn/kT) 1 / 2 and 
all times will be made dimensionless by dividing by this. 
For molecular dynamics we have to specify the masses of 
the particles. We assume that the mass of a sphere scales 
with the cube of its radius, correct if all spheres are made 
from the same substance with the same uniform density. 
The density we use is the volume fraction r\ — (ir/6)pa 3 , 
where cr 3 is the third moment of x(a). 



3 Crystallisation 

Although wc will present results only for systems which 
showed no sign of crystallising, some of our systems did 
start to crystallise. Speedy 19 has undertaken a study of 
the freezing of monodisperse hard spheres within molecular 
dynamics simulations. He finds crystallisation only above a 
volume fraction of 0.54. We simulated 1372 monodisperse 
hard spheres. They were stable for a time much longer 
than the relaxation time of the fluid at a volume fraction 
i] = 0.54 but crystallised at rj = 0.56. Thus we agree with 
the findings of Speedy. Simulations of 1372 spheres with a 
polydispersity width of w — 0.1 also crystallised but when 
the width was increased to w = 0.2 (s = 5.8%) the fluid was 
stable, it never crystallised. Simulations at polydispersity 
widths w — 0.3, 0.5 and 0.7, of both 256 and 1372 spheres, 
never showed any sign of crystallisation. They remained 
amorphous up to the highest densities simulated. Widths 
w = 0.3, 0.5 and 0.7 correspond to values of s of 8.7%, 
14.4% and 20.2%, respectively. 

Experiments on colloids find no crystallisation for poly- 
dispersities beyond a standard deviation s close to 8%. 14, 15 
Thus our finding of crystallisation when the standard de- 
viation is 3% (w = 0.1) but not when it is 6% (w = 0.2) 
is consistent with experiment. However, our simulations 
are for systems much much smaller than experiment and 
the time elapsed until nucleation occurs will depend on 
system size. Theoretical and computer simulation stud- 
ies of the equilibrium phase behaviour of polydisperse hard 
spheres 20-27 show that polydispersity destabilises the crys- 
talline phase, pushing the fluid-crystal transition to higher 
volume fractions. 

Our systems with polydispersity widths of w = 0.2, 0.3, 
0.5 and 0.7 do not crystallise even at volume fractions up 
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to 0.58 because either the fluid phase is still the equilib- 
rium phase even at these high densities or it is mctastable 
with respect to fluid-crystal coexistence but the fluid is not 
sufficiently deep into the coexistence region for the barrier 
to nucleation of the crystalline phase to be small. 3 The 
nucleation rate varies as the inverse of the product of the 
relaxation time and the exponential of the free energy bar- 
rier. 3 If the barrier is large up to densities at which the 
relaxation time is very large, say r\ = 0.57-0.58, then the 
rate will always be very small. Perhaps too small to be 
observed. 

Theoretical studies 20-22 ' 25 ' 26 suggest that for polydis- 
persity widths of around 0.3 and greater, the crystalline 
phase cannot form with the full width of the distribution. 
It can only form if it accepts only a narrow range of diame- 
ters. Neither in experiments 14 ' 15 nor here is the formation 
of a crystalline phase with a narrow range of diameters 
observed to form from a polydisperse mixture of spheres. 
It is not clear whether or not the fact that the crystalline 
phase can only form from a fraction of the distribution of 
diameters is a cause of the fact that crystallisation is not 
observed. 

We checked for segregation of the large and small 
spheres, i.e., for a tendency to phase separate into a phase 
of large spheres coexisting with a phase of small spheres. 
There was no sign of such segregation. Segregation is 
not expected at the polydispersity widths studied here in 
the fluid phase. The Boublik-Mansoori-Carnahan-Starling- 
Leland (BMCSL) theory 28,29 as generalised to polydisperse 
hard spheres by Salacuse and Stell 17 does however predict 
phase separation at much broader distributions of diam- 
eters. 30,31 This fluid-fluid separation is metastable with 
respect to phase separation plus crystallisation of the large 
spheres. 32 



4 Pressure 

The reduced pressure is plotted as a function of volume 
fraction in Fig. 1. Results from simulation are shown along 
with the predictions of the BMCSL expression for the pres- 
sure. 17 ' 28,29 The agreement between the BMCSL equation 
and simulation is excellent. The only noticeable deviation 
is that the BMCSL is a little too high at the highest den- 
sities for w = 0.3 and 0.5. 

The pressure tends to decrease as the polydispersity 
width increases. The random-close-packing density of poly- 
disperse hard spheres increases as the polydispersity width 
increases. 33,34 The pressure diverges at random-close- 
packing. Therefore, at constant r\ we move farther from the 
density where the pressure diverges as the polydispersity 
width increases and so it is unsurprising that the pressure 
decreases. The distribution used in Refs. 33,34 was differ- 
ent to the hat function used here but the trend should be 
relatively insensitive to the exact form of the distribution 
when the polydispersity is not too great. 



5 Dynamics 

As we are considering dense states the motion of a par- 
ticle is strongly restricted by its neighbours. Over some 
timescale a particle will rattle back and fore in a cage 
formed by its neighbours; then as we are considering a 
fluid not a glass, the particle will 'break out' of the cage 
and eventually its motion will be diffusive. We can assess 
whether or not particles are just rattling back and fore or 
are moving significant distances, i.e., distances compara- 
ble to their own diameter, by examining the intermediate 
scattering function, sometimes called the self incoherent in- 
termediate scattering function, F s (q,t). As indicated it is 
a function of wavevector q and time t, it is also a function 
of volume fraction and polydispersity. It is defined by 35 

1 N 

%t) = ^E cos ^- ^ (*) - r * ) < ( 2 ) 

i=i 

where rj (t) is the position vector for the ith particle at time 
t. For the fluids we consider, it is not a function of the ini- 
tial time, taken to be t = above, only of the time elapsed, 
t. For our isotropic systems it is also only a function of 
the magnitude q of the wavevector q, not of its orientation. 
To improve statistics we average over wavevectors along 
the x, y and z axes. Clearly F 8 (q,0) = 1 and at times t 
sufficiently long that the positions of the particles are no 
longer correlated with their positions at time t = 0, then 
F s (q,t) — 0. As the particles move, their positions become 
decorrelated with their positions at t = and F s decays to 
0. For a fixed q, F s (q, t) will become small when the parti- 
cles have moved a distance of about n/(2q). All our results 
are for q = 2n/a, and so F s will become small when most 
of the particles have moved a distance of around a quarter 
of their diameter. F s enables us to define a relaxation or 
correlation time for the dynamics, r. We define r by 

F s (2n/a,T) = 1/e. (3) 

Values of r for to = 0.3, 0.5 and 0.7 are shown in Fig. 2. 
As expected r increases rapidly as the volume fraction in- 
creases; by about 2 orders of magnitude between i] = 0.5 
and 0.58. It is also clear that at constant volume frac- 
tion, increasing the polydispersity decreases the relaxation 
time t. However, if the relaxation time is plotted as a 
function of the compressibility factor Z, Fig. 3, then the 
data for the 3 polydispersity widths essentially follow the 
same curve. Note that in Fig. 3 we have included data for 
monodisperse spheres (at densities below where they crys- 
tallise) and shown results down to r\ = 0.3. Spheres with 
a polydispersity width of to = 0.7 at r/ = 0.57 have almost 
the same Z and almost the same relaxation time as spheres 
with to = 0.3 at r\ = 0.56. These are the two points which 
are almost superposed in Fig. 3. To a good approxima- 
tion the relaxation time as defined by Eq. (3) is a function 
only of the compressibility factor Z = pV/(NkT). It is not 
obvious to the author why r should depend only on Z al- 
though there are theories which attempt to relate dynamic 
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quantities such as r to thermodynamic quantities. 1,3 If we 
plot t as a function of the reduced pressure p r (not shown) 
then the results for the 3 polydispersity widths do not fall 
on the same curve, in particular the results for w = 0.7 are 
above those for the narrower distributions. 

We should note that for r\ = 0.56 and above the relax- 
ation times calculated using 256 spheres are consistently 
above those calculated with 1372 spheres. For w = 0.3 
and r\ = 0.57 the relaxation times are 77 and 110 for 1372 
and 256 spheres, respectively. This suggests that the re- 
sult for 1372 spheres is 10-30% higher than that of an in- 
finite system. When t is large, we observe finite size ef- 
fects for the dynamics but not for the statics. This seems 
to be a quite general finding. 36 It suggests that there is 
some length scale associated with the dynamics but not 
the statics which grows as the relaxation time r increases. 
Recently this possibility has been extensively studied by 
Glotzer and coworkers. 12, 13,37,38 They define a purely dy- 
namical correlation length and observe that it increases as 
the dynamics slow down. Perhaps the easiest way to see 
the correlations in the dynamics is to plot out the fastest 
and slowest moving particles. We simulate the particles 
with w — 0.3 at the highest volume fraction studied for 
this polydispersity, r\ — 0.57. Fig. 4 shows the fastest 5% 
of the 1372 particles, while Fig. 5 shows the slowest 5%. 
Fastest and slowest are defined as having the largest or 
smallest displacement over a time interval of length 29.4. 
This is a little less than half the relaxation time r. Clearly, 
neither the fastest nor the slowest particles are randomly 
distributed, both show strong clustering — the dynamics 
are highly correlated. This has been observed in previous 
simulation studies. 10-13,36-38 

So over times of order r the positions of the fastest 
particles are correlated. The movement of the particles is 
very far from being independent diffusion of each particle. 
If the particles were diffusing independently then for each 
particle, the probability of finding the particle at a point 
would be a Gaussian function of the distance between the 
point and the position of the particle at t — 0. As they are 
not diffusing independently it will be not be a Gaussian. 
We can measure deviations of the distribution of particles 
from a Gaussian by using the non-Gaussian parameter in- 
troduced by Rahman 39 and defined by 

a(t)= w-EZohM ->.«»]• (4) 

5(w-'Ei 1 [r,(t)-r,(0)] a ) 

For a Gaussian distribution the ratio of the fourth moment 
to the square of the second is 5/3 and a = 0. In Fig. 6 
we have plotted a as a function of time for a number of 
volume fractions. For monodisperse particles the Maxwell- 
Boltzmann distribution of velocities enforces a — in the 
short time, ballistic, regime, and at times much longer than 
t we have diffusion and again a = 0. For polydisperse par- 
ticles a is nonzero in the ballistic regime due to different 
particles having different masses, and nonzero at very long 
times due to larger particles having smaller diffusion coef- 



ficients than smaller particles. So, in Fig. 6 a is always 
nonzero. The value of a in the ballistic limit depends only 
on the polydispersity width w, its value in the other limit 
comes from the spread in diffusion coefficients with particle 
size and so will also depend on density. Positive values of 
a come from a large fourth moment, due to more particles 
with large displacements than would be found for a Gaus- 
sian distribution with the same second moment. It is clear 
if we compare Fig. 6 with the relaxation times of Fig. 2 
that the non-Gaussian parameter a goes through a maxi- 
mum at a time comparable to but (at least at the highest 
densities) less than the relaxation time r. This maximum 
value increases rapidly as r increases. The deviations from 
a Gaussian are much larger at a volume fraction of 0.57 
than at 0.5. 

Recently, experiments have imaged colloidal suspen- 
sions of hard-sphcrc-like particles at the high volume frac- 
tions considered here. 5-7 The suspensions are of parti- 
cles with somewhat smaller polydispersities than consid- 
ered here, 8%, 5 6% 6 and 5%. 7 The dynamics of colloidal 
particles, 40 are rather different from those of free particles 
such as our hard spheres. Colloidal particles are immersed 
in fluid and so their motion is diffusive on all time scales on 
which the particles move a significant distance. In addition, 
the motion of nearby particles is coupled by hydrodynamic 
interactions through the fluid. If we compare the a values 
measured by Weeks et al. 7 with those of Fig. 6 we find that 
the simulation and experimental results are rather differ- 
ent. The experimental values of a are consistently larger 
than those in Fig. 6, Weeks et al. measure a maximum a at 
rj = 0.52 of about 1.5 and at rj = 0.56 of a little under 2.5. 
The dynamics on a time scale of the relaxation time r must 
be rather different for colloidal particles and free particles, 
although in both cases the relaxation time increases very 
rapidly and in both cases glasses can form. The results of 
Kegel and van Blaaderen 6 for a are also larger than our 
simulation results; those of Kasper et al. 5 are much larger, 
at rj = 0.56 the maximum value of a is over 5. 

Two obvious sources of the difference between molec- 
ular dynamics simulation and experiment are: i) the mo- 
tion of colloidal particles is diffusive even on time scales 
much less than t, and ii) hydrodynamic interactions be- 
tween particles, i) and ii) could be distinguished by per- 
forming both molecular dynamics and Brownian dynamics 
of the same potential. 41,42 Brownian dynamics is an ap- 
proximation to the dynamics of particles in a fluid which 
accounts for the fact that particles diffuse even over very 
small length scales but neglects hydrodynamic interactions 
between particles. 40 

Another way to look at the large spread in the amount 
individual particles are moving is to define F s ''s for each 
particle, and then use this to obtain a relaxation time for 
each particle. Averaging over wavevectors along the x, y 
and z axes we have for the ith particle 

- {cos(g[x,(C l ) - Xi(0)]) + cos( 9 [ yi (Ci) - Vi(0)])+ 
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cos(«[*i(Ci) - Zi(0)])} = 1/e, (5) 

which defines the relaxation time for the ith particle, Q. 
As always q = 2ir/a. In Fig. 7 we have plotted the proba- 
bility function P((), where P(()d( is the probability that 
a particle will have a relaxation time between ( and C + dC- 
For both r\ — 0.5 and 0.56 there is a peak at short times 
and then P decays. The decay is clearly close to expo- 
nential, which makes sense as (at not-too-short times) the 
probability of a particle leaving its cage in some short time 
interval should be independent of time; at equilibrium the 
environment of a particle or cluster of particles is on av- 
erage not changing with time. For r\ = 0.5 the peak is a 
little below r at that density but for r\ = 0.56 the peak is at 
about one tenth of r. At the higher density the probability 
distribution of ( is very broad. 

5.1 Polydispersity 

The spheres are polydisperse, different spheres have differ- 
ent diameters and so their static and dynamic properties 
will differ. We expect the smaller spheres to move larger 
distances than larger spheres in the same time. At very 
short times, in the ballistic regime, this is trivially true 
due to the smaller mass and hence higher velocity of the 
smaller spheres. We can look at motion over a time scale 
of order r by examining the sizes of the fastest and slowest 
5% of the spheres. Figs. 4 and 5 show these spheres for 
w = 0.3. We average over 8 runs with an average time 
interval close to that in Figs. 4 and 5. The average di- 
ameter of the fastest 5% (69) of the spheres, over a time 
interval of 29.2, is 0.95ct ± 0.01<7, and for the slowest 5% 
it is I.OIct ± 0.01(7. The effect is small but as wc would 
expect the fastest spheres are smaller than average. For a 
polydispersity width of 0.7, over a time interval of 5.88, the 
average diameter of the fastest 5% is 0. 79 cr ± 0.0 2 cf, and for 
the slowest 5% it is 1.08±0.02ct. A time of 5.88 is, as with 
the less polydisperse spheres, somewhat less than half r 
and so near the maximum in a. Of course the effect is now 
larger but still the slowest spheres are not much larger than 
the mean, their mean diameter would have been above 1.3 
if all the slowest spheres were also the largest. The differ- 
ence between the mean diameter of the fastest spheres and 
the overall mean diameter is twice the difference between 
the mean diameter of the slowest spheres and the overall 
mean diameter. 

We have defined an intermediate scattering function, 
F s , for all spheres regardless of diameter in Eq. (2). We 
can also define an intermediate scattering function for sub- 
sets of the particles with diameters in some range <7 m j„ 
to a m ax, F s (q,t;a min ,a max )- The definition is completely 
analogous to that of the total F s (q,t) in Eq. (2), 

— — cos(q.[r 4 (t)-r,(0)]), 

(6) 



where the sum is over all spheres with diameters in the 
range a m i„ to cr max , and N mm is the number of spheres 
in this range. As before q — 2it /a and we average over 
wavevectors along the x, y and z axes. We will study in- 
termediate scattering functions for the smallest, Fg S \q,t), 
and largest spheres, Fs l \q,t), defined by 

F$'\q,t) - F s (g,t;<f(l-«V2),<f(l-u;/2 + 0.05)) 
F^(q,t) = F s (q,t;a(l + w/2-0.05),a(l + w/2)). 

(7) 

These are for the particles within 0.05(7 of the minimum 
or maximum diameter. For w = 0.7 this corresponds to 
spheres with diameters in the ranges 0.65 to 0.7(7 and 1.3 
to 1.35(7. We also define relaxation times t^ s ' and t"> from 
Fs^ (q, t) and Fjf 1 (q, t), respectively, using the analogues of 
Eq. (3). t( s ) and t''' for w = 0.7 are plotted as filled trian- 
gles in Fig. 2. The ratio r' 1 '/^ 8 ' increases with increasing 
density but slower than r does. At r\ = 0.57 the ratio is 
approximately 10. 

A feature of the dynamics when r is large is that they 
are heterogeneous, some particles travel much farther than 
others in a time r. In a fluid of monodisperse particles 
these heterogeneities are dynamic, a particle may be fast 
at one time but later on may be slow. For polydisperse 
particles in addition to this source of heterogeneity there is 
the spread in diameters which means that some particles, 
the smaller ones, are on average faster than others over all 
time scales. In Fig. 6 we can compare the a parameter 
for w = 0.3 and 0.7 at the same volume fraction r\ = 0.57. 
a is larger for the more polydisperse spheres; this is the 
case despite the fact that the relaxation time for the more 
polydisperse particles is smaller. The system with w = 0.7 
and r\ = 0.57 is at almost the same compressibility fac- 
tor as the one with w = 0.3 and i] = 0.56 so we see that 
at constant compressibility factor a increases sharply with 
polydispersity. 

For weakly polydisperse spheres the relaxation times for 
all sizes of sphere increase together and so if we are at high 
enough density then r will be very large and on a timescale 
much less than r all spheres will be localised, none or at 
least very few will have moved more than a small fraction 
of their diameter. However, if spheres with a wide range of 
diameters are present then the smallest particles will have a 
relaxation time scale which is much less than r. Then over 
some time scale much less than r but much larger than the 
relaxation time scale for the smallest particles not all the 
spheres will be localised. The smallest spheres will be dif- 
fusing. This has been observed in binary mixtures of hard 
spheres in which the smaller spheres are much smaller than 
the larger spheres, by Jackson et al 43 The same effect is 
seen in crystals of binary mixtures. 44 This is an extreme 
example of polydispersity making the dynamics even more 
heterogeneous than they are when the particles are all the 
same size. Note that polydispersity and the heterogeneous 
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nature of the dynamics when r is large have the same effect, 
they produce a wide spread in speeds/relaxation times of 
the individual particles. Even for monodisperse particles, 
for times much less than r the heterogeneous dynamics 
mean that some particles have left their cages. However, 
if we have some particles much smaller than others then 
the small ones may have a relaxation time not much larger 
and a diffusion constant 43 not much smaller than at lower 
densities. 

6 Discussion 

We have studied the dynamics of dense polydispcrsc hard 
spheres using molecular dynamics. The relaxation time r at 
a given volume fraction was found to depend on the polydis- 
persity, it decreased as the polydispersity increased. So, the 
glass transition is pushed to higher volume fractions as the 
spheres become more polydisperse. This is consistent with 
the increase in the density of random-close-packing with in- 
creasing polydispersity. 33 ' 34 Although we have found this 
using molecular dynamics, we expect that experiments on 
polydisperse colloids would show that as polydispersity in- 
creased the kinetic glass transition observed in experiment 
would move to higher volume fractions. 

The relaxation time r is, to a good approximation, a 
function only of the compressibility factor Z, changing the 
polydispersity at fixed Z has no effect on r. This is despite 
the fact that a characteristic of the dynamics, such as a, 
changes a great deal at fixed Z as w increases; compare 
the w = 0.3, 77 = 0.56 and w — 0.7, rj — 0.57 curves in 
Fig. 6, which have almost the same Z but very different 
a's. Apparently, increasing the polydispersity at constant 
volume fraction reduces the relaxation time by reducing Z. 
The virial equation, 45 
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relates Z to the forces between the particles. () denotes an 
ensemble average. In Eq. (8) and fy are the vector be- 
tween the centres of the ith and jth particles and the force 
on the ith particle due to the jth particle, respectively. 
The sum is over all pairs of particles. So, Z is a constant, 
1, plus a term proportional to the sum over the product 
of the interparticle forces and the interparticle separations. 
Why the relaxation time should be a function only of this 
product is not clear to the author. 

Although Z can be expressed in terms of forces it is also 
of course a thermodynamic quantity. At constant tempera- 
ture Z varies as p/p so r being a function of Z is equivalent 
to it being a function of the ratio of the pressure to the 
number density. As the energy of hard spheres is purely 
kinetic it is independent of volume, thus the pressure is 
simply the volume derivative of the entropy times the tem- 
perature. So, the relaxation time is a function only of Z 
and so is a function only of the ratio of volume derivative of 
the entropy to the number density. But it is not a simple 



function of it; r does not vary as the exponential of the 
Z, it increases more rapidly. The pressure shows no sign 
of a discontinuity in slope so we conclude that there is no 
phase transition in the density range studied here. It is an 
open question whether or not there is a phase transition 
to an 'ideal' glass phase at higher densities, as assumed by 
Speedy, 4 see also Refs. 1,3. 

We have characterised the dynamics using the non- 
Gaussian parameter a, Fig. 6, and the distribution of 
relaxation times P(C), Fig. 7. The results for a have 
been compared to experimental results for colloidal suspen- 
sions. 5 The dynamics are heterogeneous 10-13 by which 
we mean that there is a broad distribution of speeds and 
of relaxation times of the particles. This is true not only 
at times of order the relaxation time r but also, at higher 
density and hence r, for times at least an order of magni- 
tude larger, see Fig. 7. Indeed given the very wide spread 
of individual relaxation times £ it is clear that the overall 
relaxation time r is inadequate to characterise the time de- 
pendence of the dynamics. In particular, relaxation is far 
from complete even at times much longer than r. Com- 
parison of our non-Gaussian parameter with that obtained 
in experiments on colloidal hard spheres showed significant 
differences. One possible explanation of this is that hydro- 
dynamic interactions between colloidal particles are acting 
to make the relaxation more cooperative, i.e., acting to in- 
crease the clustering shown in Fig. 4. It is known that 
hydrodynamic interactions tend to favour collective over 
relative motion. 40 
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Figure 1: The reduced pressure of polydisperse hard spheres. The solid, dashed and dotted curves are the predictions 
of the BMCSL theory for polydispersity widths of w — 0.3, 0.5 and 0.7, respectively. The circles, squares and triangles 
are the results of simulation, for polydispersity widths of w = 0.3, 0.5 and 0.7, respectively. 




Figure 2: The reduced relaxation time r as a function of volume fraction 77. The open circles, squares and triangles are 
the results of simulation, for polydispersity widths of w = 0.3, 0.5 and 0.7, respectively. The filled triangles are and 
r^, the relaxation times of the the smallest and largest spheres, at a width of w = 0.7. The smallest spheres have the 
smaller r at all densities. 
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Figure 3: The reduced relaxation time r as a function of the reduced pressure Z. The circles, squares and triangles arc 
the results of simulation, for polydispersity widths of w = 0.3, 0.5 and 0.7, respectively. The +'s are for monodisperse 
spheres. For w = 0.3 and 0.7, and for monodisperse spheres we have extended the calculations down to i] = 0.3. 
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Figure 4: The fastest 5% of our particles, 69 out 1372, are drawn as spheres of diameter a (all the spheres are drawn the 
same size). Fastest is defined as having traveled the largest distance within a time interval of 29.4. The volume fraction 
tj = 0.57 and w = 0.3. 
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Figure 5: The slowest 5% of our particles, 69 out 1372, are drawn as spheres of diameter a. Slowest is defined as having 
traveled the smallest distance within a time interval of 29.4. The volume fraction rj = 0.57 and w — 0.3. 




Figure 6: The non-Gaussian parameter a as a function of time. The circles and squares are for a polydispersity widths 
w = 0.3 and 0.7, respectively. The three curves for w = 0.3 are from bottom to top, for r\ = 0.52, 0.56 and 0.57. The 
squares are for 77 = 0.57. 
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Figure 7: The probability distribution function for relaxation times of individual particles. The crosses and squares are 
for rj = 0.5 and 0.56, respectively; w — 0.3. 
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